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i ABSTRACT 

J>^' Using high resolution N-body simulations of stellar disks embedded in cosmologically 

I motivated dark matter halos, we study the evolution of bars and the transfer of an- 

. gular momentum between halos and bars. We find that dynamical friction results in 

' some transfer of angular momentum to the halo, but the effect is much smaller than 

! previously found in low resolution simulations and is incompatible with early analyt- 
ical estimates. After 5 Gyrs of evolution the stellar component loses only 5% - 7% of 

' its initial angular momentum. 

, Mass and force resolutions are crucial for the modeling of bar dynamics. In low 

■ resolution (300 - 500 pc) simulations we find that the bar slows down and angular 
I momentum is lost relatively fast. In simulations with millions of particles reaching a 

j. . resolution of 20-40 pc, the pattern speed may not change over billions of years. Our 

' high resolution models produce bars which are fast rotators, where the ratio of the 

, corotation radius to the bar major semi-axis lies in the range TZ — 1.2 — 1.7, marginally 

■ compatible with observational results. In contrast to many previous simulations, we 
' find that bars are relatively short. As in many observed cases, the bar major semi-axis 

"q^, is close to the exponential length of the disk. 
I ' The transfer of angular momentum between inner and outer parts of the disk plays 

O I a very important role in the secular evolution of the disk and the bar. The bar forma- 

^ , tion increases the exponential length of the disk by a factor of 1.2 -1.5. The transfer 

^ ■ substantially increases the stellar mass in the center of the galaxy and decreases the 

\ dark matter-to-baryons ratio. As the result, the central 2 kpc region is always strongly 

^ . dominated by the baryonic component. At intermediate (3-10 kpc) scales the disk 

' is sub-dominant. These models demonstrate that the efficiency of angular momentum 

rS I transfer to the dark matter has been greatly overestimated. More realistic models 

■ produce bar structure in striking agreement with observational results. 

Key words: galaxies: kinematics and dynamics — galaxies: evolution. 



1 INTRODUCTION 

Making predictions for the structure of individual galax- 
ies in the framework of cosmological models is a dif- 
ficult task, that until recently was mostly dominated 
by modeling of dark matter dominated systems such as 
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well on those tests predicting too fast rotation in the 
central region. Excessive amount of dark matter satellites 
is an additional complication for the cosmological models 
jKlvpin et alJll999l:lMoore et alJll999D . 

Another layer of problems is on slightly larger scales: 
luminous parts of normal galaxies and masses in the range 



10^^M„ . The number of issues is quite large. 



An excessive amount of dark matter in th e bulge may 
result in too few micro-lensing events (e.g.. IZhao fc Mao 



19961: iBinnev fc Evand 1200 J) , but see also jKlvpin et al 



20021) . Mass modeling of the Milky Way galaxy on dif- 



ferent scal«s_js_a^_conmlicate and quite uncertain issue 
(e.g.. IPehnen fc Binnevlll998fl . For halos predicted by cos- 
mologic al models the conclusions change from very pes- 
simistic ("n avarro fc Steinmet dl2000l :l Hernandez et al .1 1200 II ) 
to ne utral ( Eke et al, 2001) to more optimistic (.Klypin et alJ 
l2002i). 



One of the contentious issues is the pattern speed of 
bars. Using analytical arguments and a combination of N - 
body simulations with r otatin g solid bars IWeinberd il985l) : 
iHernauist fc Weinberd ^1992^ argued that a non-rotating 
dark matter will exert so much tidal friction that a bar 
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should lose its angular momentum in few rotation periods. 
For a rigid rotating bar with mass comparable to the mass 
of the spheroid (includi n g dar k matter) inside the bar radius 
iHernauist fc Weinberd lll992fl predict that the bar "sheds 
all its angular momentum typically in less than ~ 10^ yr". 
If true, that would make impossible for bars to exist in 
galaxies with dark matter halos. There are no doubts that 
the dynamical friction operates when a bar rotates inside a 
dark matter halo. Yet, the rate and the amount of the an- 
gular momentum lost by the disk is still a matter of debate. 
Fully self-consistent N-body si mulations with live bars , 
disks, and dark matter halos fCombes fc S anderj Il98lt 
ISellw o od fc Athanassoula 1986; Dcbattist a fc SellwoodI 
ri998l."l20oJ lAthanassoula fc Misiriotij I2OO2I) show that 
bar s slow down far less than what has b e en pr edicted 
by IWeinberd ^8^; IHernauist fc Weinberd lll992ll . Bars 
in simulations rotate for billions of years, defying the 
theoretical expectations. Y et, bars in numerical models 
are o bserved to slow down. iDebattista fc SellwoodI il998l 
|200[|) find that in models, in which initially the dark 
matter constitutes about the same amount of mass as the 
stellar disk in the central part of the galaxy, the baryonic 
component (disk -I- bar) loses about 40 percent of it angular 
momentum over ~ 10 Gyrs. The bar pattern speed declines 
by a factor of five after the bar forms. This produces a bar 
which rotates more than twice slower than the observed 
bars. The speed of bar rotation is often characterized by 
the ratio TZ of radius of corotation to the large semi-axis of 
the bar. 

Pattern speed s for real gala xies can be es timated us- 
ing the method of iTremaine fc W einberg lll984) . For a few 
galaxies - most of them are SO' s - the ratio was measured 
and is in the range TZ 1.1-1.7 iMerrifield fc Kuiikenll995l: 
Gerssen et alJll999l: lAguerri et alJl2003l) . 3ars in models of 
Debattista fc SellwoodI lll998l . bOOOiT liave 7^ = 2.0 - 2.6, 
which contradicts the observed values. 

A less k nown but not less important problem of bars is 
their length. ICombes fc ElmegreerJ (^93) compare lengths 
of bars in N-body simulations with those observed in galax- 
ies. Even though a good agreement with observations was 
found, the models either did not include a live halo or in- 
cluded only a bulge component and no halo at all. In more 
recent simulations with live halos (dark matter dominated 
or not) the bar semi-major axes ar e 2-4 times longer than 
the i n itial d isk exponential length ( Debattist a fc SellwoodI 
Il998l. I2OOOI: lAthanassoula fc Misiriotis 2002). These long 
bars contradict observations because real galaxies typically 
have bar lengths 0.5-1.5 of the present time exponential 
disk scale iElmegreen fc Elmegreeiilll985|). For example, the 
bar in our Galaxy is 3-3.5 kp c long ( Blit z fc SpergeJll99lt 
IZhadll996l: lFreudenreichlll99^ : ICerhard 2003) which is close 
to the disk scale length of 2 . 5-3.5 kpc jPehnen fc BinnevI 
ll998l : lHammerslev et al . '"l999';"Gerhard"200:^. In the case of 
the models studied by jDebattista fc Sollwood (1998, 2000) 
the bar extended over the whole disk. That would cor- 
respond to 15 kpc, if we scale the models to match our 
Galaxy. Scaling the dark matter dominated model MH of 
lAthanas soula fc Misirioti^ (|20o3) to our Galaxy, the radius 
of the bar would be 14 kpc (about 4 times the initial expo- 
nential length). 

The arguments - the slowing down of bars and the 
lengths of bars - against significant amount of the dark 



matter in the central parts of normal galaxies have prob- 
lems on their own. The arguments are based on numeri- 
cal simulations, which suffer from two types of problems: 
numerical effects and initial conditions. Numerical resolu- 
tion is very limited in most simulations. It is often em- 
phasized that models should have a large number of par- 
ticles to sup press the two-body scattering and to reduce th e 
noise (e.g., iweinberd Il99i IDebattista fc SellwoodI I2OO0I1 . 
Yet, the force resolution is another important effect, which 
is often overlooked. For example, the formal (one cell) res- 
olution in simulations of IDebattista fc SellwoodI ll2000h was 
1/5 of the disk scale length. The scale hight was 1/2 cell, 
which means that vertical structure was not resolved. If 
we scale the disk scale length to 3.5 kpc of our Galaxy, 
then we get the formal resolution of 700 pc. To make the 
things w orse, the real resolution i n the Particle Mesh code 
used by IDebattista fc Sellwood' f2000|) is not the formal 
(cell) resolu tion, but about twice the formal resolution. For 
IDebattista fc SellwoodI 1119981) and IDebattista fc SellwoodI 
i200(]() simulations this corresponds to 1.4 kpc, which is 
clearly not enough to mim ic th e disk of our Galaxy. The 
resolution in simulations of IFuxI (Il997ll was 1.8 kpc at the 
solar distance. 

There are number of reasons why high resolution (both 
mass and force) is needed for simulations of the bar for- 
mation. (1) The central region of galactic models. Propa- 
gation of waves through the center and swing amplifica - 
tion are prime suspects fo r bar instabilitv iToonir^ il98lll : 
iBinnev fc Tremaind (I^S^)- Simulations must have sufficient 
resolution to allow correct treatment of these waves. High 
resolution at the center is also required because during the 
bar format ion the density at the center inc reases very dra- 
matically (lAthanassoula fc Misirioti^ l2002l . e.g.). This in- 
fall of the stellar component to the center should be ac- 
curately tracked: the distribution of mass in the central re- 
gion may affect the structure of the bar ( Norin an et al.-199^ : 
iBerentzen et alJll998tlShen fc Sellwoodli20oll . (2) The ver- 
tical structure of the disk. Disks of galaxies are very thin - 
only a couple hundred kpc. Resolving the vertical structure 
of disks is a challenge for numerical simulations: even opti- 
mistically the resolution should be not worse than 100 kpc. 
Tracking the vertical evolution of the disk is even more dif- 
ficult than one naively expects. At some stages of the evo- 
lution t he disk develops waves which o s cillate in vertical di- 
rection (' Fridman fc Pohachenkol Il984l: iMerritt fc SellwoodI 
Instabilities related with these waves (e.g. the fire- 
hose inst ability) are often blamed for sudden thickening o f 
the disk dRaha et al.ll99ll:lAthanassoula fc Misiriotisl2002ll . 
In this case the code should be able to maintain both the 
vertical structure and to treat sufficiently accurately the up- 
ward and downward displacements of the thin disk. The ver- 
tical oscillations combined with random radial stellar veloc- 
ities produce coupling of the vertical and the radial direc- 
tions. As the result, high resolution in the vertical direction 
must be accompanied with the same resolution in the ra- 
dial direction. (3) Non linear coupling of waves. The role of 
interactions between different modes is poorly un derstood 
(IContopoulojll98ll : iFridman fc Pohachenkol Il984^ . Theory 
operates with line ar perturbations and typic ally ignores the 
nonlinearities (e.g. iNelson fc Trema ine"l999Y This does not 
mean that the nonlinearities are not important: we simply 
do not know how to handle them analytically. There are no 
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doubts that these effects exist: bars have very large density 
contrasts. As such, they are nonhnear features. Even more 
difficult problems are related with the reaction of the dark 
matter to the small-scale effects in the disk. The dark mat- 
ter reacts to what happens in the disk sometimes even in a 
resonant fashion llAthanassoulalEnO^Tt . In turn, disk reacts 
to the disturbances in the halo. This non-linear interaction 
is clearly important for the long-term evolution of the bar 
iAth anassoulal l2003l) . Numerical simulations are the only 
available tools to study the effects. Coupling of long and 
short waves (e.g. the short vertical frequencies with long 
m — 2 modes) require very high resolution. It should be 
noted that imposing a high force resolution without ade- 
quate mass resolution can be a disaster: simulations would 
amplify small-scale shot noise potentially leading to incor- 
rect predictions. 

Initial setup is questionable in many cases when it 
comes to the properties of the dark matter halo. If one 
wants to test predictions of cosmological models, those pre- 
dictions should be used for setting initial conditions of nu- 
merical models. For example, halos predicted by cosmologi- 
cal models are very large and should e xtend to 200-300 kp c 
for galaxies of the size of our Galaxy jKlvpin et al.ll2002l) . 
Halos i n simulatio ns are much smaller than that. For ex- 
ample. |Fux ( 1997) us ed halo trunc ated at 38 kpc. Mod- 
els of iDebattist a fc SellwoodI (l200dl were larger - 40 kpc, 
but still too small as compared to the virial radius of the 
halo of our Galaxy. The density profile in the outer part 
of a halo (radii larger 20 kp c for the Milky Way galaxy) 
should decline as pdm oc iNavarro et aljll997f) . This is 
very different from the profile s used in m any simulations of 
bar formation. For example, IRd] ^1997^ used exponential 
profile Pdm cx: expf— r /9.1 kpcTi lAthE massoula fc Misiriotij 
^2002) used density profile p oc exp(— r^/(35 kp c)^)/[r^ 4- 
(1.75 k pc)^] truncated at 52.5 kpc. IDebattista fc SellwoodI 
(|2000j) used polytrops n = 3, which have little relation with 
the expected dark matter profile. 

The large size and mass of cosmological dark matter ha- 
los affect the motion of halo particles even in the central re- 
gion. Because more massive halos produce deeper potential 
wells and have larger escape velocities, dark matter particles 
move faster in the central region. In turn, larger velocities of 
particles make t he halo more resistan t to interactions with 
the stellar disk. lAthanassoulal (|200^ argues that "hotter" 
halos absorb less angular and result in slower decline of the 
bar pattern speed. 

Both problems of the existing numerical simulations - 
the resolution and the initial conditions - motivate us to 
simulate with high resolution models, which use more real- 
istic initial conditions. We use the traditional approach to 
form a bar: bars are formed by dynamical instability i n ini- 
tially featureless equilibrium disk rotating inside halos fHohi 
Il971: Ostrikor & Peebl es 1973: Miller 1978). In section|2|we 
describe initial conditions for our models and give details 
of the code used to run our simulations. The evolution of 
models is discussed in section 2] In section |5| we present 
results on the angular momentum of different components. 
The pattern speed and the structure of bars are presented 
in section |5| We discuss our results and compare them with 
the previous results in section |7| A brief summary of our 
conclusions is given in section |H| 



2 MODELS 

2.1 Initial conditions: densities and velocities 

We study models, which initially have only two components: 
an exponential disk and a dark matter halo. No initial bulge 
or bar are used. Initia l cond itions for the systems are gener- 
ated using iHernauistl il993l) method. We use the following 
approximation for the density of the stellar disk in cylindri- 
cal coordinates: 



Pd{R,z) 



sech (z/zo), 



(1) 



where Md is the mass of the disk, Rd is the exponential 
length, and zq is the scale height. The later was assumed to 
be constant through the disk. The disk is truncated at five 
exponential lengths and at three scale heights: R < 5Rd, 
z < 320. Disk mass was slightly corrected to include the 
effects of truncation. The vertical velocity dispersion is 
related to the surface stellar density E and the disk height 
zo- 



gI{R) = ■kGzoY.{R), 



(2) 



where G is the gravitational constant. The radial velocity 
dispersion or is also assumed to be directly related to the 
surface density: 



(Tr{R) = Ae 



(3) 



The normalization constant A in eq. Q is fixed in such a 
way that at some reference radius Ry-d the radial random 
velocities are Q times the critical value needed to stabilize 
a differentially rotating disk against local perturbations: 



(^b{Rt. 



, 3.36GS(i^i., 

K(_Rrcf) 



(4) 



where Ac(f?ref) is the epicycle frequency at the reference ra- 
dius. The reference radius is chosen to be -Rrof = 5 kpc for 
our models. Results are very insensitive to the particular 
choice of the reference radius because the stability parame- 
ter Q does not change much for radii from Rd to few Rd- We 
also consider a family of models where Q is fixed along the 
disk. In this case eq.l|lj gives the radial dispersion at every 
point. The fast increase of the epicycle frequency k in the 
central region overcomes the exponential growth of density 
making the disk colder in the centre. 

The tangential (</>) component of the rotational velocity 
and its dispersion are found using the asymmetric drift and 
the epicyclic approximations: 



2 



2R 
Wd ^ 



_ 2 K 



(5) 



(6) 



Here the circular velocity Vc is found as the quadrature sum 
of the halo and the disk contributions. For the disk compo- 
nent we use the thin-disk approximation. 

We assume that the da rk matter density pr ofile is de- 
scribed by the NFW profile (iNavarro et alJll997f) : 



Pdn 



x(l + a;)2' 



r/rs 



ln(l + C) - 



C 



1 + C 



(7) 



C=— , (8) 
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where Mvir and C are the virial mass and the concentration 
of the halo. For given virial mass the virial radius of the halo 
is found assuming a flat cosmological model with matter 
density parameter flo = 0.3 and the Hubble constant Ho = 
70 kms"^Mpc"\ 

Knowing the mass distribution of the system M{r), we 
find the radial velocity dispersion of the dark matter: 

1 r GMjr) , 

(^r.dm = / Pdm 5 dr.. (9) 

The other two components of the velocity dispersion are 
equal to cr^^^. In other words, the velocity distribution is 
isotropic, which is a good appr oximation for the c entral parts 
of the dark matter halos fe.g.. lColin et al.ll200ol) . Eq.® ig- 
nores non-spherical deviations of the gravitational accelera- 
tion for the dark matter. At each radius we find the escape 
velocity and remove particles moving faster than the escape 
velocity. 

The dark matter halo is truncated at the virial radius, 
which for our models is at 250 kpc - 300 kpc radius. Dur- 
ing the evolution of the system, the truncation results in a 
gradual smearing of the outer boundary of the halo. Some 
particles move to radii larger than the virial radius. They 
are not replaced by other particles, which, in true equilib- 
rium, would come from the region outside the truncation 
radius. This creates a rare-faction wave moving inside the 
halo. The wave never reached the central part of the halo 
and was always outside the central 100 kpc region. 

The distributions eqs. H1I9^ are used to produce initial 
coordinates and velocities of particles. The space is divided 
into bins. We estimate the expected number of particles in- 
side each bin. We than randomly place specified number of 
particles inside the bin. Spherical shells are used to initi- 
ate dark matter particles. For the disk we use cylindrical 
shells along radius and equal-spaced bins along the axis of 
rotation. Velocities are picked from an appropriate Gaussian 
distribution truncated at the escape velocity. The number of 
particles was typically 10-20 in each bin. The resulting dis- 
tribution of particles was tested against expected analytical 
expressions. 

The disks are realized with particles of an equal mass. 
The dark matter halo is composed of particles of variable 
mass with small particles placed in the central region and 
larger particles at larger distances. We double the particle 
mass as we move further and further from the centre. In total 
we use 5 mass species with the mass range of 32. The region 
covered by small mass particles is ~ 40 kpc - much larger 
than the size of the disk. In order to reduce the two-body 
scattering each dark matter particle in the central region 
has the same mass as a disk particle. 

This procedure is designed to reduce the number of par- 
ticles and, at the same time, to allow us to cover a very large 
volume. For example, in model Ai we have 3.55 million par- 
ticles inside ~ 300 kpc, of which 2.4 million are small parti- 
cles. We would need 9.5 million particles, if particles of equal 
mass were used. In the course of evolution some of massive 
particles come closer to the centre, but their number was 
always very small. 



2.2 Choice of parameters 

The approximation eq. © formally is valid only for pure 
dark matter halos without baryons. Contraction of baryons 
in the process of formation of the disk should significantly 
affect the dark matter. We try to mimic the effects of the 
contraction by changing the halo concentration and by tun- 
ing parameters of the halo and the disk in such a way that 
the distribution of the dark matter rough ly approximates 
more realistic models. iKlvpin et al.l l|2002) presented such 
models for the Milky Way and for the M31 galaxies. The 
models used realistic cosmological dark matter halos and 
satisfied numerous observational constraints. To some de- 
gree we mimic two favored models in that paper: models Ai 
and Bi, which are described in detail in Klvpin et al. (200^. 

The first model (Ai) was designed assuming that the 
angular momentum of the dark matter was preserved during 
baryonic collapse. The conservation of angular momentum 
leads to a specific prediction that in the central ~ 5 kpc the 
density of the dark matter closely follows the density of the 
baryonic component. The second model (-Bi) includes effects 
of an exchange of angular momentum between the baryons 
and the dark matter. The dark matter gains angular mo- 
mentum through the dynamical friction during early stages 
of disk and bulge formation, when non-axisymmetric pertur- 
bations are expec ted to be large . Some what similar scenario 
was proposed bv lEl-Zant et al] i200ll) . Model Bi also had 
significant amount of the dark matter in the central region, 
but less than in model Ai. For example, the crossing point, 
at which the contribution of the dark matter to the circular 
velocity is equal to that of the baryons, is shifted from 5 kpc 
in model Ai to 12 kpc in model Bi. 

In order to reproduce the effects of the contraction and 
the exchange of the angular momentum, we increase the con- 
centration of the dark matter halos. Ou r models have C = 1 5 
as compared with C = 12 for models of lKlvpin et al.H2002l) . 
We also adjust parameters of our models (virial mass, ex- 
ponential length and mass of the di s k) to roughly mimic 
important properties of iKlvpin et alJ (l2002h models in the 
central 20 kpc. For example, our models {Ai and A2) have 
a sub- maximum disk (contributions of the disk and the halo 
to the circular velocity are approximately equal) at radii 
R < 5 kpc with halo dominating at larger radii. The mass 
of the disk was 4.28 x lO^^'M ^ , which is c l ose to the mass 
of disk+bulge in model Ai of iKlvpin et alj J2002l) . 

The only difference between our two models Ai and A2 
is in the velocity dispersions of disk particles. For model 
Al we use eq.l^l with the normalization provided by eq.Q 
at the reference radius. In the case of model A2 the eq.l^J 
is used at all radii (not just at one radius) and the eq.lO 
is not used at all. In other words, model A2 has the same 
Q parameter at all radii. The main difference between the 
models is that model A2 is "colder" in the central region. In 
the central ~ 1 kpc region model A2 has significantly smaller 
radial and azimuthal velocity dispersions. Both models have 
the same circular velocity and the same vertical velocity 
dispersion. 

Our next model {B) has a smaller halo and a more mas- 
sive disk (6 X IO^^Mq ) . The disk mass is equal to the sum of 
the di sk and the bulge masses in model Bi of lKlvpin et alJ 
i2002h . We have a third model called C. This model has a 
disk mass similar to that in models A, but both radial and 
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Figure 3. Initial (dashed curves) and final (full curves) rms ve- 
locities of the stellar component for the model C. The left panel 
shows the radial and the right panel shows the vertical (z) com- 
ponents of the velocity. Initially the rms velocities are exponential 
with the ratio of cTrad/i^vert ~ 2.5. Formation of the bar "heats" 
the vertical component more than the radial one. Outside the bar 
(R > 5 kpc) most of the change in the rms velocities occurs in the 
first 1 Gyr when spiral arms form and dissipate. Notice that in 
this region the radial and azimutal components of velocities are 
preferentially heated by the spiral arms. 



Figure 1. Initial and final circular velocities for model Ai (left 
panels) and for model B (right panels). Dashed curves show 
the contribution of the stellar component GM{r)/r. The dot- 
dashed curves are for the dark matter. Final models show the to- 
tal stellar contribution: disk and bar contributions are combined. 
Initially model Ai has a sub-maximum disk with the density of 
the disk being close to the density of the dark matter inside the 
central 3 kpc region. Model B is more dominated by the disk 
component: the crossing point of the disk and dark matter is at 
9 kpc. 




2 4 6 8 10 12 14 2 4 6 B 10 12 14 
R (kpe) R (kpc) 



Figure 2. Initial (left panel) and final (right panel) circular veloc- 
ities for model C. Dashed cu rves show the contribution of the stel- 
lar component -y/ GM {r)/r. The full and the dot-dashed curves 
are for the total circular velocity and for the contribution of the 
dark matter correspondingly. 



vertical scale lengths are shorter. The halo has the same 
mass as halo in model Bbut concentration is considerable 
larger. 

Parameters of our models are presented in Tabled Top 
p anels in Fi gure ^ shows the initial circular velocities Vc = 
^J'GMir)Jr of the disk and the halo for the models. 

There are important differences between our models and 
the models described in lKlvpin et alJ i2002h . The later are 
realistic models, which, in addition to the exponential disk 
and the halo, have also a nucleus and a triaxial bar. The goal 
of our paper is to study how a bar develops in an unstable 



Table 1. Initial parameters of models 



Parameter 




B 


C 


Disk Mass (10i°Mg ) 


4.28 


6.0 


4.8 


Total Mass (lO^^Mg ) 


2.04 


1.0 


1.0 


Disk exponential length (kpc) 


3.5 


3.0 


2.9 


Disk exponential height (kpc) 


0.25 


0.25 


0.14 


Stability parameter Q 


1.2 


1.3 


1.2 


Halo concentration after 


15 


15 


19 


baryonic contraction C 








Total number of particles (10^) 


35.5 


7.8 


97.7 


Number of disk particles (10^) 


2 


1.1 


12.9 


Particle mass (IO^Mq ) 


2.14 


5.45 


0.37 


Maximum resolution (pc) 


22 


43.6 


100. 



disk. Thus, our models initially do not have either a bar or 
a bulge: all the stellar mass is in the exponential disk. The 
bar and the bulge are expected to develop f rom the disk. 
Initially, our models mimic to some degree iKlvpin et al.l 
(2002) models. For example, the exponential disk length is 
3 — 3.5 kpc. Yet the subsequent evolution of the models was 
quite significant. Thus, the final parameters a rather dif- 
ferent from what was assumed at the beginning. Our final 
models do not provide a fit to the Milky Way galaxy. Nev- 
ertheless, the gross properties (e.g., the maximum circular 
velocity and the stellar surface density at 8.5 kpc) are not far 
from what is observed in our Galaxy. This is why we scale 
different physical quantities to astronomical units (kpc, , 
and so on). 



3 NUMERICAL SIMULATIONS 
3.1 Code 

We use the Ada ptive- Refinement-Tree (ART) A^-body code 
iKravtsov et al. 1997; Kravtsov 1999) to run the numerical 
simulations analyzed in this paper. The code starts with a 
uniform grid, which covers the whole computational box. 
This grid defines the lowest (zeroth) level of resolution of 
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the simulation. The standard Particles-Mesh algorithm is 
used to compute the density and gravitational potential on 
the zeroth-level mesh with periodical boundary conditions. 
The code then reaches high force resolution by refining all 
high density regions using an automated refinement algo- 
rithm. The refinements are recursive. A refined region can 
also be refined. Each subsequent refinement level has half 
of the previous level's cell size. This creates a hierarchy of 
refinement meshes of different resolution, size, and geometry 
covering regions of interest. Because each individual cubic 
cell can be refined, the shape of the refinement mesh can be 
arbitrary and effectively match the geometry of the region 
of interest. This algorithm is well suited for simulations of a 
selected region within a large computational box, as in the 
simulations presented below. 

The criterion for refinement is the local density of par- 
ticles. If the number of particles in a mesh cell (as estimated 
by the Cloud-In-Cell method) exceeds the level rithrosh, the 
cell is split ("refined") into 8 cells of the next refinement 
level. The refinement threshold depends on the refinement 
level. The threshold for cell refinement was low on the zeroth 
level: nthrosh(O) — 2. Thus, every zeroth-level cell containing 
two or more particles was refined. The threshold was higher 
on deeper levels of refinement: nthrcsh = 3 and nthrcsh ~ 4 
for the first level and higher levels, respectively. 

During the integration, spatial refinement is accompa- 
nied by temporal refinement. Namely, each level of refine- 
ment is integrated with its own time step, which decreases 
by factor two with each refinement. This variable time step- 
ping is very important for accuracy of the results. As the 
force resolution increases, more steps are needed to integrate 
the trajectories accurately. 

The ART code has the ability to handle particles of dif- 
ferent masses. We use multiple masses to increase the mass 
(and correspondingly the force) resolution inside a region 
centered around the simulated galaxy. 



3.2 Parameters of simulations 

The simulations presented here are run using f28'' zeroth- 
level grid in a computational box of 2.86 Mpc for models 
^and Sand 1.4 Mpc for the model C. The models are placed 
in the center of the grid and far from the periodic images. 
Because the size of the models are small as compared with 
the size of the box, the effects of the periodical images are 
very small and can be neglected. For example, at a distance 
of 100 kpc from the center of the model ^the contribution of 
periodic images is less that 7 x 10~^ of the main galaxy force. 
Even at 250 kpc, which is close to the virial radius for the 
models, the tidal force is only 6 x 10~^ of the force from the 
central image. Bar formation and evolution are practically 
unaffected by the periodic boundary conditions. 

Model Ai is simulated for 8.2 x lO' yrs with 550 zero- 
level time steps. The code reaches 10 refinement levels and 
has about 9 million cells. The number of time steps at the 
highest resolution level is about 560,000, which gives the 
time step of 1.5 x 10* yrs. With this time step it takes 3500 
steps to make one orbital period at the 3 kpc radius. The 
formal spatial resolution - a cell size at the 10th level - is 
22 pc. At about twice that value the force is close to the 
Newton's law. Central 6 kpc region was resolved with cells 
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Figure 4. Test of the code: The growth of the power spectrum 
of perturbations in the standard cosmological model with 512'^ 
particles. The dashed curves show the theoretical power spectra at 
different moments. The full curves show results of the simulation. 
The wiggles of the simulated spectrum at small wavenumbers k 
are due to the small number of long- wave harmonics. They should 
be preserved (shifted upward). From bottom to top the curves 
are for redshifts z = 48, 18.8, 9. At z = 9 the fluctuations in the 
simulation have slightly larger amplitude at large wavenumber 
due to nonlinear effects. The code accurately tracks the growth 
of about a million of independent harmonics when the amplitude 
of fluctuations increases many times. 

of the size not more than 100 pc. Model A2 has similar 
parameters, but it was simulated for 7 Gyrs. 

Model B is simulated with a smaller time step of 1.02 x 
10"* yrs (at the 9th level) for 4.5 x 10^ yrs. It reached the 
highest resolution of 43.6 pc, but because of fewer particles, 
the volume with high resolution is significantly smaller than 
in models A: there are only 3.9 million cells in model B. The 
number of steps at the highest 9th level is 4.4 x 10^. 

Model Chas the largest number of particles: almost 10 
million of which more than a million particles belong to the 
disk. The resolution of this simulation was forced to be not 
better than 100 pc. Because of large number of particles, the 
whole disk (up to 14 kpc in radius) is resolved with 100 pc 
cells. The time step at that resolution was 1.2 x 10^ yrs. The 
simulation box was 1.4 Mpc. One of the motivations for the 
model C is to significantly reduce discreteness effects. For 
example, our estimate of the two-body scattering time is 
more than 10^ Gyrs for this model (see sectiorj^jlj. 

It is interesting to compare our numerical simu- 
lations with those presented in ^ebattista & Sollwoo^ 
Jl998l I2OO0I) and with more recent simulations of 
lAthanassoula fc Misirioti^ i2002h . It is straightforward to 
compare the resolutio n of o u r sim ulations with that in 
iDebattista fc SellwoodI Jl998l |2000D because we use the 
same Particle-Mesh algorithm. The formal (one cell) reso- 
lution in their simulations was 1/5 of the disk scale length. 
If we scale the length to 3.5 kpc, than we get the formal 
resolution of 700 pc, which is 30 times larger than in our 
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models. lAthanassoula fc Misiriotij J2002ll used the TREE 
code with the Plummer softening of 0.0625 of the disk ex- 
ponential length. Again, scaling it to 3.5 kpc, we get 220 pc. 
The force of gravity comes close to the Newton's law at ap- 
proximately three Plummer softening lengths or at 660 pc. 
In our models A this scale is 40-80 pc in the centre and is 
200 pc for the entire 6 kpc region. 

3.3 Tests of the code 

Extensive tests of the code and compa risons with othe r nu- 
merical A''-body codes can be found in'Kr avtsovl il99!j) and 
fKnebe et al. 120001 .Kravtsov (1999 ) (Figure 6) shows the 
accuracy of the gravitational force. The code was tested 
against known analytical solutions: (1) the growth of small 
perturbations in the expanding universe; (2) the spherical 
accretion ("the Bertschinger solution"); (3) the collapse of 
a plane wave (the Zeldovich solution). In all tests the code 
performed extremely well. For example, in the case of the 
spherical accretion the density at the center should increase 
by 5 orders of magnitude over tested period of time. The 
code reproduced that density increase as well as the whole 
radial density profile. Here we present results of two addi- 
tional tests, which are more relevant for the bar models: 
the growth rates of perturbations and the stability of a self- 
gravitating equilibrium system. 

Figure 01 shows the evolution of the power spectrum of 
perturbations in the standard cosmological model. We use 
512"^ particles to set initial conditions at very high redshift 
2; = 50. The dashed curves in the Figure 01 show the theo- 
retical power spectra at different moments. The full curves 
show results of the simulation. At the initial moment the de- 
viations of the simulated spectrum from the analytical one 
at small wave numbers k are due to a small number of long- 
wave harmonics. This is normal. The code should preserve 
those deviations at later times, which it clearly does. This 
test shows that the code is able to accurately follow the evo- 
lution of numerous waves over extended period of time dur- 
ing which the waves increase their amplitude many times. 
For a code this is a very difficult test. It probes every block 
of the code: gravity solver, motion of particles, and den- 
sity assignment. There are numerous independent harmonics 
involved (millions in our case). Each wave must grow inde- 
pendently on others in a well defined way. Because the long 
waves have much larger amplitude, any erroneous small cou- 
pling with short scales would ruin the power spectrum. The 
code was able to accurately track the growth of all the waves 
during a period when amplitude of the waves increases many 
times. 

Tests also were performed to check the accuracy and 
stability of bounded orbits of particles. For example, we set 
two particles in a circular motion around each other and 
and run the code with the two particles for hundreds of 
orbits. We do not find any drift of radius of the orbit. Even 
a small systematic loss or gain of angular momentum or any 
adverse effect when a particle moves from one refinement 
level to another would be clearly detected by the test. The 
code was tuned not to have any. 

We present results of a test of this kind: a long-term 
stability of an equilibrium self-gravitating configuration. We 
use 200,000 particles to set an equilibrium halo with the 
NFW density profile. Three mass species of particles are 



Table 2. Parameters for the test Halo Simulation 



Parameter Model 


Hi 


Halo Mass (Mg ) 


1 X 10^2 


Halo concentration 


15 


Total number of particles 


2 X 10^ 


Total number of mass species 


3 


Formal maximum resolution (pc) 


174 



used, with the second specie being twice as massive as the 
first one and the third one being four times as massive as 
the first specie. The halo was truncated at the virial radius 
of 250 kpc. Parameters of the halo simulation are presented 
in Table 121 The simulation is then run for 5 Gyrs. There are 
different goals for this test. Effects of the two-body scatter- 
ing are discussed in the next section. Here we present the 
evolution of the density profile. Figure (51 shows the initial 
and final profiles of the halo. The only observed change in 
the profile is the decline of the density close to the virial 
radius. This is caused by the truncation of the configuration 
at the virial radius. A rarefaction wave gradually moves in 
to inner parts of the halo. Yet the effect of the wave is very 
small even at large radii. For example, the density at 100 kpc 
decline only by 10 percent. The central part (r < 60 kpc) of 
the halo does not show any measurable change of the den- 
sity. Note that the crossing time for the central 2 kpc region 
is only lO^yrs for typical velocity 200 kms~^. Thus, the cen- 
tral part was run for 500 crossing times. The NFW profile 
is a difficult test because of the large densities and low ve- 
locity dispersion in the central region. This cold and dense 
center potentially is unstable: any numerical errors have a 
tendency to heat up the center and to reduce its density. 
The code passes this tough test. 

3.4 Two-Body Scattering 

Two-body scattering should be very small for systems, which 
we try to model. Given that our halo and disk are made of 
particles, every precaution should be made to insure that 
discreteness does not affect the evolution of the systems. 
On small scales this is done by imposing a low limit on the 
number of particles inside a resolution element. The force of 
gravity becomes Newtonian at about 2 cells. The code is de- 
signed to have at least 4 particles in a cells for models A and 
B and about 20 for the model C. Thus, inside one resolution 
element - a sphere with a radius of 2 cells - there are 130 and 
670 particles for models A,B and C correspondingly. This 
large number of particles means that the force acting on a 
particle can never be dominated by its few closest neighbors. 
In other words, the code totally eliminates the most dam- 
aging close collisions. They cannot possibly happen in the 
code we use. 

Large distance encounters are possible and really hap- 
pen in our models, but their effect is small. We use the sim- 
ulation Hi to investigate the role of the two-body scatter- 
ing in our immerical experiments and, specifically, to study 
whether or not angular momentum is affected by this pro- 
cess. The angular momentum for each individual particle 
should be preserved, if the model is truly collision- less. The 
change of the orbital angular momentum of any given i— th 
particle, ALi = |Li — Li,init|, could be attributed to numer- 
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Figure 5. Test of the code: stability of an equilibrium NFW halo 
modeled with 200,000 particles for 5 Gyrs. The dashed and the 
full curves show the initial and the final profiles. The central part 
of the halo is potentially susceptible to different numerical effects 
and may be affected by the two-body scattering. The profile inside 
100 kpc shows no evolution. The slight decline of density outside 
100 kpc is due to initial truncation of the halo at virial radius of 
250 kpc. 



leal effects or to the two-body scattering or to a combina- 
tion of both. If it is dominated by the two-body scattering, 
the random walk approximation is a reasonable model for 
the change of orbital angular momentum (Chandrasekhar 
1942). In this case the rms value of the change of the orbital 
angular momentum -y/EALf/iV should grow as the square 
root of time. The relaxation time associated with the two- 
body scattering can be defined as the time when the rms 
change in the orbital angular momentum is comp arable to 
the rms value of the orbital angular momentum -y/ELf/iV. 
This could be expressed as: 



EL? 



(10) 



The right hand size of ea. llUl is a linear function of 
t with the slope equal to the inverse of relaxation time 
(TRciax). Figure |S| shows the behavior of EAL^/E(I/i) over 
5 Gyrs. The linear behavior shown in the plot supports the 
random walk approximation expressed in ea. llUII . The slope 
corresponds to a relaxation time of 3.2 x 10^ Gyrs. As a 
test we can compare our result with theoretical estimates. 
Chandrasekhar (1942, equation 2.380) gives an analytical 
expression for Tkciax due to two-body scattering: 

((73d/20)3 



1.12 X lO" 



,A 



Gyrs, 



(11) 



nimim2 log^p . 

where mi and 7712 are the mass of the colliding particles 
in solar masses, 0-313 is the three dimensional rms velocity 
of the scattered particles in kms"'^, ni the number density 
of the scattering particles in pc~^ and In A is the Coulomb 
logarithm. An estimate of Tadax based on a Foker-Planck 
orbital average, which also uses velocity deviations instead 



of energies, predicts a coefficient for eg. 1111 that differs only 
by a factor of 0.6. We neglect this small correction. In the 
case of scattering between particles of the first mass specie 
(the dominant term) we should use in eg. lllll the mass and 
the number density of small particles. For the Coulomb loga- 
rithm one should take InA = ln(bmax/&min), where 6niax and 
femin are the maximum and minimum values for the collision 
impact parameter. A reasonable approximation for these val- 
ues is: bniax ~ Rvir ~ 252 kpc and femin = force resolution — 
0.174 kpc. This gives InA = 7.3. We calculate the value of 
the relaxation time using ea. dllt . If we take r = 40 kpc ( 25 
% mass radius), n = 1.5 x 10"^°pc-^, mi = 2 x 10^ Mg , 
(T3D = 171 km/s, we obtain TRciax = 3.6 x 10'^ Gyrs, which 
is consistent with the estimate obtained from our halo sim- 
ulation. 

A disk-halo system modeled with particles potentially 
may experience a large two-body scattering if there is a large 
difference in masses between particles. This happened in old 
simulations, when massive halos where modeled with rela- 
tively few particles to reduce the computational time. Our 
main simulations A- C use 5 mass species, each mass specie 
being twice massive as the previous one. Thus, the mass 
ratio of the most massive to the lightest particle is large 
and equal 16. In the initial configuration the central 40 kpc 
of the halo and of the disk is sampled only with the par- 
ticles of the same (small) mass. This setup is designed to 
minimize the artificial scattering. However, few particles of 
the 5th (most massive) specie eventually reach the central 
region. How does this affect the two-body scattering relax- 
ation time estimate? An estimate of the relaxation time Tbig 
for particles of the 5-th specie giving energy to ("heating") 
the particles of the 1-st mass specie (small particles) can 
be calculated using ea. ljll^ by replacing the terms nimim2 
with nfcmfcms, where indexes b and s refer to big and small 
particles correspondingly. Thus, the relaxation time of the 
1st specie due to scattering with big particles scales as 



±b ~ -[Rcla 

ribmb 



(12) 



Even at late moments of evolution the number density of 
massive particles in the central 1 kpc is 4 orders of magni- 
tude smaller than the number density of the lightest par- 
ticles. Using those numbers we obtain that Ti, is 625 times 
longer than the relaxation time for the first specie TReiax- 
There is an additional effect, which increases this estimate 
even more. Eq. 1121 does not take into account the fact that 
large particles move a factor of two faster when they move 
through the central region as compared with the small par- 
ticles. This happens because large particles have large ener- 
gies: initially they are placed at large radii. In other words 
the two-body relaxation due collisions with big particles is 
negligible in our simulations. 

As an additional test for the effects of the two-body 
scattering we study how the halo density profile evolves 
during the simulation: two-body scattering should lead to 
a gradual decline in the central density for the NFW halo. 
Figure |^ shows the evolution of the density profile of the 
halo simulation after 5 Gyrs. The central region of the halo 
has a higher density and a shorter dynamical time. How- 
ever we do not find any change in the density profile even 
after 5 Gyrs. This confirms our estimates of the relaxation 
time-scale. 
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Figure 6. The test of the code for preserving the angular mo- 
mentum of individual particles. The plot shows the time evolution 
of the average deviation of the angular momentum AL; = |Li — 
Li initi in the Hi simulation with 2 X 10^ particles (crosses). The 
line is the AL oc \/T fit for the points after 0.5 Gyro. The slope 
implies the two-body relaxation time Tj^eiax = 3.2 X 10^ Gyrs, in 
agreement with the Chandrasekhar formula evaluated at 30 kpc. 

The simulation Hi is used just as a test for code sta- 
bility and for estimates of the two-body scattering. Our 
main simulations have many more particles, and, thus the 
scattering should be even smaller for them. Using eq. lllll 
and assuming virial equilibrium we obtain a scaling rela- 
tion: Taoiax oc iV/ln7V (Binney & Tremaine 1987). This 
relation allows us to scale-up the estimates for the relax- 
ation time based on simulation Hi to our simulations with 
a larger number of particles. For our disk-halo simulations 
with 3.5 X 10® and 10^ particles we get TReiax = 4.5 x 10'' 
Gyrs and 1.2 x 10^ Gyrs respectively. We can safely ignore 
any effects of scattering in our disk-halo simulations. 



3.5 Finding the bar 

We follow iDebattista fc SellwoodI J2000D and 
lAthanassoula fc^^igirl^nn200^ when finding bars in our 
simulations. We start with finding the centre of the system. 
This c an be done in different ways. Debattista fc Sellwood 
i2000l) minimized the sum of distances of all particles 
relative to a centre. We used a more simple algorithm: 
the centre of mass of stellar particles. We then bin stellar 
particles using cylindrical shells spaced equally in logarithm 
of the distance. For each bin we find the Fourier components 
of second harmonic of the angular distribution of particles: 

JV JV 

a2 = j^J2 ««"(2</'0, ''2 = ^1] cos(20O, (13) 
1=1 1=1 

where A'' is the number of stellar particles in the bin and (ja 
is the angle of the i— th particles. We then find the phase 
and the amplitude of the second harmonic A2 = (a^ + b2)/2. 



Except for the very early stages of evolution, when the 
bar was just emerging, it is easy to identify the bar. In 
the central region the amplitude of the second harmonic A2 
changes gradually from bin to bin and its phase is almost 
constant. At larger distances (typically larger than 5 kpc) 
the amplitude declines and the phase shows large variations. 
In the very centre the shot noise complicates the results. 
Because of these considerations, we use bins with radii in 
the range 0.5-3 kpc to find the amplitude and the phase 
of the bars. We find the bar amplitude and its phase as the 
amplitude-weighted means at different radii: (A2) — yv^f), 
(0) = (0^2) / (^2)- When finding the average amplitude and 
the phase, we exclude bins with A2 less than 1/2 of the max- 
imum amplitude found in the range of the radii. Once the 
phase of the bar is found at different moments of time, we 
estimate the frequency of bar rotation (the pattern speed) 
by numerical differentiation: flp = d(j>/dt. 

Finding the radius of the bar is a more difficult 
problem because bars do not have clear boundaries. Dis- 
cussion of different methods of finding the bar major 
ax is can be found in [Debattis ta 8z SellwoodI i2000l) and 
in lAthanassoula fc Misirioti^ ~ il200^ . We use two different 
methods, each being sensible and producing reasonable, but 
slightly different results. We find the surface density along 
and perpendicular to the bar. We then define the radius of 
the bar as the radius at which the surface density along 
the bar is twice of the surface density at the same dis- 
tance perpendicular to the bar. In the second prescription 
the bar radius is defined as the radius at which the con- 
tours of the surface density have axial ratio of two. We also 
tried other prescriptions, w hich are to some degree rem i- 
niscent to what was used bv IDebattista fc SellwoodI J200CI) . 
We find the major semi-axis as the radius at which either 
the phase of the bar starts to deviate by more than 10 de- 
grees or as the radius at which the amplitude falls below 
1/2 of its maximum value. We find that typically the pre- 
scriptions based on the surface density give the same re- 
sults as those based on the amplitude and the phase of the 
bar. For example, for model A2 &t t — 6.2 Gyrs the bar 
radu are: J?bar = 4.3 kpc (A2,max/2), = 5.5 kpc (A0 = 10°), 

— 4.2 kpc (Smajor (^major ) — 2Sniinor (^major ) ) , — 5.3 kpC 
(Smajor ('^major ) — ^minor (^major /2) ) . Here Smajor (^) and 

Sminor('") are the stellar surface densities along the major 
and the minor axes respectively. Nevertheless, we find that 
the surface densities give more stable results because occa- 
sionally spiral waves in the outer disk align with the bar and 
the amplitude/phase conditions produce spurious results. 

We present the total range of the bar sizes produced by 
the two methods based on the surface density profiles. The 
estimates vary by 20-30 percent. We would like to empha- 
size that every estimate was reasonable: the estimates in the 
middle of the range are no better than the extremes. 



4 EVOLUTION 

Figure Q shows the evolution of the stellar component in 
model Ai. Because the system is relatively hot and the dark 
matter is significant in the central region, the bar devel- 
ops slowly. After 1 Gyr (top left panel) the bar has not 
yet appeared. Nevertheless, extended spiral arms are clearly 
present in the disk. At that time the system has already 
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Figure 7. Distribution of the stellar component at different stages of evolution for model A-^. Starting from the left bottom panel in the 
clock-wise direction the panels correspond to the initial moment; 1 Gyr; 2 Gyr; 3.3 Gyr. The disk rotates counter-clockwise. To avoid 
crowding we show only every fifth particle. At the distance 5 kpc the disk made 20 orbital periods during 3.3 Gyr. 



made many rotations: at the distance of 5 kpc one orbital 
period is equal to T « 1.6 x 10® yr. The bar develops af- 
ter 2 Gyrs (top right panel). At that moment its radius is 
« 4 — 5 kpc. The spiral arms can still be seen, but they start 
to get weaker after that. At the last shown moment (3.3 Gyr) 
only fragments of spiral arms can be detected. The bar is 
still very strong, but its amplitude is visibly smaller than at 
the peak at 1.8 Gyr. The bar gradually grows with time. At 
later moments (5 8 Gyrs) the major semi-axis of the bar is 
« 6 kpc. The evolution of model B is qualitatively the same, 



but it's initial evolution is more violent and is faster than in 
the case of model A]_. 

The evolution of model A2 remarkably differs from that 
of model Ai in spite of the fact that their initial parameters 
arc almost the same. The only difference between the models 
is in the initial random stellar velocities: the central part of 
the disk is colder in model A2. The initial evolution is faster 
in model A2. After 1 Gyr the bar is already in place. The 
bar is shorter (1 - 1.5 kpc) and it rotates almost twice faster 
then the bar in model Ai. The bar very slowly increases its 
length. This slow evolution proceeds for 2 Gyrs. At 3 Gyrs 
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Figure 8. Evolution of fraction of disk mass inside given radius. 
The top panel is for models Ai (full curves) and A2. The bottom 
panel shows model B. The total mass inside 5 kpc practically 
does not change with time. Most of the mass exchange happens 
inside the central 1-2 kpc. The disk mass inside the central 
1 kpc region increased by a factor 4 by the end of evolution. The 
time-scale of evolution is about twice shorter in the case of the 
more disk-dominated model B and the colder model ^2- 



from the beginning the bar suddenly starts to grow in length 
and its amplitude also increased significantly. 

Formation of bars is accompanied by significant redis- 
tribution of the stellar mass. Even a visual comparison of 
the bottom panels in Figure |7| shows that the central sur- 
face density increases. Figure |H| shows the evolution of the 
fraction of stellar mass inside given radius. The total mass 
inside 5 kpc practically does not change with time. All the 
mass exchange happens inside the central 1-2 kpc and it is 
very violent. For example, the stellar mass inside the cen- 
tral 1 kpc region increases by a factor 4 (5.5) after 5 Gyrs 
(8 Gyrs) of evolution. Note that the time-scale of the ini- 
tial evolution is about twice shorter in the case of the more 
disk-dominated model B and for the colder model A2. 

Bottom panels in Figure |3| present more detailed infor- 
mation on the change in the densities. In the central 2 kpc 
the ratio of the stellar density to the dark matter density 
have increased by a factor of two by the end of the evolution. 
As the result, even model Ai became significantly baryon- 
dominated in the centre, the same effect was already ob- 
served in simulations by lAthanassoula & Misiriotis (2002). 
The Figure also shows that the central growth of the density 
goes at the expense of the stellar density in the middle part 
( 2-8 kpc ). At those radii the final baryon-to-dark matter 
ratio was a bit surprising: it went down producing a dark- 
matter dominated disk. This effect is more pronounced in 
model B, which was initially more baryon-dominated than 
models A. 

The dramatic changes of the disk component hardly 
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Figure 9. The baryon-to-dark matter ratios (bottom panels) and 
the ratio of the rotational velocity to the velocity dispersion (top 
panels) for model Ai (left panels) and model B (right panels). 
The dash curves show the initial models. The final models are 
shown with full curves. The central growth of the density goes at 
the expense of the stellar density in the middle 2 — 8 kpc region. 
At those radii the final baryon-to-dark matter ratio went down 
producing dark-matter dominated disk. The top panels show that 
the heating of the disk preferentially happens in the outer regions 
outside of the bar. Spiral waves may cause the heating. 

affect the dark matter: it stays almost unchanged. There is 
a 20 (40) per cent increase of the dark matter density in 
the central 2 kpc region after 5 (8) Gyrs of evolution . The 
rest of the halo is practica l ly un affected in agreement with 
lAthanassoula fc Misiriotid (|20o3). 

Because of the violent processes in the central region, 
one would naively expect to find a very strong heating of 
the disk in the central part and much less in the outer parts. 
The ratio of the rotational velocity V,^ to the velocity dis- 
persion characterizes the "hotness" of disk. The top plots 
in the Figure]^ show just an opposite trend: the heating is 
relatively stronger in the outer part of the disk. A possi- 
ble explanation for this is the heating by spiral waves that 
appeared in the outer parts of the disk. Large amplitude 
spiral waves produce non-circular motions and increase the 
azimuthal dispersion. After the waves die out, they leave 
behind a hotter disk. 



5 ANGULAR MOMENTUM 

Formation of a bar is accompanied and even may be driven 
by the exchange of angular momentum between different 
components of the system. The exchange of angular mo- 
mentum between the quickly rotating stellar component and 
the non-rotating dark matter is of special interest because 
is often used as one of the arguments against the presence 
of a significant amount of the dark matter in the central 
parts of galaxies. Numerically it is a challenge to simulate 
the two-component system: any artificial numerical coupling 
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will result in a transfer of some angular momentum from the 
stellar to the dark matter component. For an axisymmetric 
system no transfer should happen. Once the bar or any other 
non-axisymmetric perturbation is formed, the dark matter 
should react to it by generating a wake behind the pertur- 
bation, which must result in a transfer of the angular mo- 
mentum to the dark matter. The exchange of the angular 
momentum is not limited by the disk-dark matter interac- 
tion. The stellar component also experiences the exchange 
between different parts of the disk, which we find to be much 
more important for the evolution and for the structure of 
bars. 

Figure \TU[ presents the evolution of the z-component of 
the angular momentum of the stellar component (disk + 
bar; top panels), the pattern speed of bars (middle panels), 
and the bar amplitudes (bottom panels) for all the mod- 
els. Figure [TTI shows details of the evolution of the angular 
momentum for Ai and B models. The total stellar angular 
momentum clearly is declining with time, but the rate of the 
decline is very small. After 5 Gyrs the change was ~ 5% for 
all the models. The disk of the model Ai lost 13 per cent of 
the angular momentum after 8 Gyrs of evolution. This level 
of the angular momentum loss is at odds with the results 
of iDebattista fc SellwoodI (j^OOO), who find about 40% de- 
cline. A large fraction of the angular momentum loss in the 
model Ai happened during a short period (5-6 Gyrs after 
the start), when the bar was increasing its amplitude. When 
the growth abruptly stopped at t ~ 6 Gyrs, the rate of the 
decline also slowed down dramatically. 

It was interesting to find that the specific angular mo- 
mentum Lz at different radii in Figure [TTI shows very little 
evolution. Only at the very centre the angular momentum 
declines slightly. At larger radii the specific angular momen- 
tum is practically constant. At first sight this indicates that 
the distribution of the angular momentum does not evolve. 
This is not true. Figure |H| shows that the stellar mass in 
the central region increases quite substantially. Because the 
specific angular momentum (r) at a given radius does not 
change, this mass increase implies that the stellar particles 
that accumulate in the central region, actually lose their an- 
gular momentum. This is clarified by Figure [T^ that shows 
significant changes in the distribution of the angular mo- 
mentum of the stellar particles in the course of the evolution. 
The changes are especially noticeable at low : a large peak 
forms. The peak at very low (almost zero) angular momen- 
tum is made of the particles located in the central region. 
The particles initially had intermediate angular momentum 
and lost some of it during the bar formation. As the result 
of this migration of particles the distribution of at later 
moments has a deep at intermediate angular momenta. The 
number of particles at large also increases, which is more 
clearly seen in the case of model B. 

Evolution of the Lz distribution is not gradual. Most 
of the changes occur before and during the bar formation. 
Model B clearly shows that once the bar has formed (after 
1 Gyr of evolution) the distribution experiences very little 
evolution. This indicates that the bar itself does not lead 
to a significant exchange of the angular momentum between 
the particles. Most of the evolution is done not by the bar, 
but by the processes that produce the bar. Nevertheless, the 
evolution does not stop after the bar forms. It proceeds, but 
with a much lower rate. 
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Figure 11. Evolution of the specific angular momentum of the 
disk for the Ai (top panel) and the B (bottom panel) models. The 
top curves in each panel show that the total angular momentum 
of disk particles declines very little — few per cent after 4.5 Gyrs 
of evolution. The other curves present an average specific angular 
momentum for disk particles inside spherical shells indicated in 
the plot. Only the very central region (radius less than 1 kpc) 
exhibits some decrease (a factor of 1.5) of angular momentum. 
Outside the very centre the specific angular momentum does not 
change. 



Because of the large moment of inertia of the dark mat- 
ter, it is difficult to detect small changes in the distribution 
of the angular momentum of the dark matter particles. The 
largest effect was found in the central region. By the end of 
the evolution the dark matter in the central 2 kpc region 
rotated with velocity 3-7 kms~^. This rotation has very 
little dynamical effect. Figure ITHl presents an example of the 
change of the angular momentum of the dark matter par- 
ticles. It shows the distribution of the z— component of the 
specific angular momentum of 10,000 particles, which ini- 
tially were in a shell centered at 3 kpc from the centre. The 
change in the angular momentum is clearly detected, but it 
is very small. 

To summarize, we clearly find indications of the dynam- 
ical friction between the dark matter and the stellar compo- 
nent. Nevertheless, the amount of the angular momentum 
lost by the stellar particles is very small: about 5% during 
5 Gyrs of evolution. There is much larger exchange of the 
angular momentum between stellar particles. Most of this 
exchange happens during the formation of the bar. Once the 
bar is in place, the distribution of the angular momentum 
changes very little. 
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Figure 10. Evolution of the bar amplitude (bottom panels), the pattern speed (middle panels), and disk angular momentum for the 
models Ai (left column of panels), A2 (middle column of panels) and B (right column). In all cases the angular momentum declines due 
to the dynamical friction with the dark matter. The decline is very small - only few percent for ~ 5 Gyrs of evolution. Note the complex 
behavior of bar amplitude and its correlation with the pattern speed. The pattern speed remains almost constant during periods, when 
the bar amplitude does not evolve (e.g., t = 2 — 5 Gyrs for model Aior t = 5 — 7 Gyrs for model A2). Periods of fast increase of the bar 
amplitude correlate with the decline of the pattern speed (e.g.,t = 3 — 5 Gyrs for model ^2). 



6 STRUCTURE OF BARS: THE PATTERN 
SPEED AND THE SURFACE DENSITY 

Figure [TC)1 shows the evolution of the pattern speed and the 
amplitude of the bars in our simulations. Models Ai and B 
evolve qualitatively similar. After different periods of time 
models Ai and Bform a bar which rotates with almost con- 
stant speed for about 5 Gyrs. The amplitude of the bar in 



the model B experienced a noticeable decline (a factor of 
two) at 2 Gyrs. This change in the amplitude was not ac- 
companied by any other obvious changes. For example, the 
pattern speed did not change. 

The bar pattern speed is slightly larger in model 
B (« 33kms~^/kpc) as compared with model Aj (~ 
27kms~^/ kpc). What is remarkable about the pattern 
speeds is that they do not change much. For exam- 
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Figure 12. Distribution of the z-componcnt of the angular mo- 
mentum of the stellar particles for models Ai (top panel) and B 
(bottom panel) at different moments of time. Initial distribution 
is shown by dot-dashed curves. Particles with low angular momen- 
tum are preferentially at small radii. Large angular momentum 
is coming from particles at large distances. Dashed curves are for 
1-1.5 Gyrs after the beginning of the evolution. The final distri- 
bution (full curves, 4.5-5 Gyrs) is qualitatively the same for both 
models. The peak at small angular momenta corresponds to the 
bar. The number of particles with intermediate angular momenta 
substantially decreases during the evolution. There is an excess 
of particles with very large angular momentum. The changes in 
the distribution are much stronger in the case of model B, which 
has less dark matter and has a more massive disk. 



Figure 13. Distribution of z-component of the angular momen- 
tum of dark matter particles at different moments. Dashed and 
full curves are for initial and advanced moments (3.3 Gyrs for 
Ai and 4.5 Gyrs for B) of evolution. We present results only for 
10000 particles, which initially were in a spherical shell centered 
at 3 kpc. Particles at other radii showed even smaller effects. 
The changes in the angular momentum are clearly seen in both 
models: there are more particles with high angular momentum. 
This indicates that the dynamical friction transfers some angular 
momentum to the dark matter, but the effect is extremely weak. 

a change even in the outer part: the exponential scale-length 
increases quite substantially - by a factor of 1.5. 



pie, in model Ai the bar changed its frequency from « 
27kms~^/kpc when it formed to « 20kms~^/kpc at t = 
8.5 Gyrs. 

The evolution of model A2 is very different. For the 
first 3 Gyrs the bar was very quickly rotating with the pat- 
tern speed of ~ 50 km s~^ /kpc - almost twice larger than in 
model Ai. At later moments the pattern speed was declining. 
It went down to ~ 25kms^^/kpc by the end of evolution. 

In order to find how fast is the bar rotation, we show in 
the Figure im the bar pattern speeds, the rotation, and the 
circular velocity curves for the models. The thick horizontal 
bars in the plots indicate the length of the bars. The corota- 
tion radius can be identified as the radius of intersection of 
the bar pattern speed with the rotation curve. In all three 
models the bars rotate fast with the ratio of the corotation 
radius to the bar length equal to 7?. « 1.2 — 1.7. 

Figure \TE\ shows the surface density profiles of the stel- 
lar components in the models. The surface density profiles 
seems to indicate a development of normal galactic compo- 
nents: an exponential disk and a (exponential) bulge. Note 
that even after significant evolution in the central part, the 
outer part of the disk (radii larger than 5 kpc) is still well 
approximated by a simple exponential profile. Still, there is 



7 DISCUSSION 

Bars in our mo dels are different from what was found 
in simulations of iDebattista fc SellwoodI (120001 . In two of 
our dark matter-dominated models bars do not show a 
strong tendency to slo w dow n. This contradicts results of 
IDebattista fc SellwoodI (■2000) , who always find that the 
pattern speed of bars substantially declines with time. We 
tried to find an explanation for the disagreement. Numer- 
ical effects are the obvi ous first suspe cts. T he force res- 
olution in the IDebattis ta & ScUwood, i200(]|) models was 
more than an order of magnitude lower than in our sim- 
ulation. In an effort to reprod u ce at least some conclu- 
sions of IDebattista fc SellwoodI J2000t) . we re-run model 
Ai with the resolution comparable with that in models of 
^obattista & Sollwood (2000J. The evolution was so much 
affected by the lack of the resolution that after running the 
model for more than 2 Gyrs and not finding even a hint 
of a bar, we abandoned the simulation. We then re-started 
the simulation with decreased "thermal" velocities in the 
dis k. The Q para meter wa s set a t the value Q = 0.05 used 
bv iDebcitti sta & S ellwoodI i200(]|) . In this simulation a bar 
forms. Figure lT^ shows that the pattern speed decreases very 
substantially over a relatively short period of time - more or 
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Figure 14. Different frequencies in models Ai, A2, and B. The 
full curves show the frequency of rotation of the disk f2rot = 
V^jR. The dashed curves present the frequency of rotation for 

a cold disk ficirc = 

The difference between the two 
frequencies are due to the asymmetric drift. The dotted horizontal 
lines show the frequency of rotation of the bars. The extension 
of the bars in the models is indicated by the thick lines. The 
distance of a bar from the rotation curves is a measure of how 
fast is the bar. The bar in the model B is clearly a very fast bar, 
which extends almost up to the corotation. The bars in models A 
rotate slower, but they are still fast bars. 



Figure 15. Surface density profiles of the stellar components for 
models A\, A2, and B. The full curves present the azimuthally av- 
eraged profiles. Dashed curves show double exponential fits with 
parameters presented in the plots. The dotted curves show the 
surface density along the major axes of the bars. The surface 
densities along the minor bar axis are shown with the dot-dashed 
curves. All the curves are obtained by averaging over three time 
moments covering about 0.5 Gyr. For comparison a large dot 
shows the surface density of the disk (stars -I- gas) for our Galaxy 
at the solar position. 



less in line with what was found bv IPebattista fc SellwoodI 
l|2000l) . Thus, we find the same behavior of bars - declining 
pattern speed - if we substantially reduce the resolution and 
reduce the random velocities of disk particles. The angular 
momentum of the stellar particles also was declining faster 
than in the high resolution run. This exercise clearly demon- 
strates that the force resolution and the stellar velocities are 
important factors. In our simulations the low force resolu- 
tion produced erroneous results that the bar pattern speed 
declines with time. The "coldness" of the disk may also have 
affected the pattern speed. This has also been reported by 
Athanassoula (2003). 

Length is another important property of bars. In 
our simulations bars have radii comparable to both ini- 
tial and final exponential disk length. This is quite dif- 
ferent from what, for example, Athanassoula & Misiriotis 
l|2002h find. The bar length in their models is more than 
three times longer than the initial disk scale-length. Af- 
ter 8.5 Gyrs of evolution (t ~ 610 in units used by 
lAthanassoula fc Misirio tis ^2002*)) the bar length in our Ai 
model is 6 — 6.5 kpc as c ompared with 10 — 14 kpc for model 
MH in lAthanassoula fc Misiriotis (2002). Bars are als o very 
long in the models of iDebattista fc SellwoodI i2000l) : they 
span the entire disk. Even that a comparison with disk scale- 
length after bar formation could be more favored for those 
models it is interesting to test how numerical resolution may 
affect the bar length. We rerun model A2 with a formal res- 
olution of 350 pc w hich is 1.6 times the formal resolution of 
lAthanassoula fc M isiriotis (2002). Figure DTI shows the dis- 



tribution of stellar particles of the simulations for model A2 
at different moments of time. The system was clearly af- 
fected by the resolution. At 2 Gyrs the bar is barely visible 
in the high resolution run, but it is much longer and stronger 
in the low resolution run. The differences are not as large 
at later moments, but even at later moments the differences 
are substantial. This example shows that the bar length is 
quite sensitive to the effects of the resolution. S hort bars 
are al so formed in high resolution simulations of ISellwoodI 
J2002I) . 

We may only speculate why a better force resolution 
leads to shorter and weaker bars. Orbital resonances are tra- 
ditionally blamed for the growth of bars. Fast rotating bars 
trap particles. This in turn increases the strength of the bar. 
It gets longer and more massive trapping even more trajec- 
tories. Low force and high mass resolutions make the system 
more resonant by producing a gravitational potential with 
too shallow gradients. Motion of particles in such potentials 
can be more susceptible to resonances. Thus, we get longer 
bars. Increased force resolution results in steeper gradients of 
the density and the potential. As the result, particle trajec- 
tories do not stay as much in resonances, which decreases the 
rate of the growth of the bar. A central conce ntration makes 
bars weaker and can even destroy them (e.g.. lNorman et alJ 
|l996?). During first stages of bar formation the central con- 
centration increases very dramatically (see Figures 8 and 9). 
A low force resolution artificially makes the central concen- 
tration shallower. This may result in longer and stronger 
bars, which are known to rotate slower. The comparison of 
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Figure 16. Evolution of the pattern speed and the bar amplitude 
in a simulation with initially cold disk with Q = 0.05. Other 
parameters of the simulation are the same as for model Ai. The 
simulation was done with a very low resolution of 0.5 kpc to 
mimic results of Debatista & Sellwood (2000). In this case the 
pattern speed quickly decreases with time giving an impression 
that the bar is slowed down by the dynamical friction with the 
dark matter. This "slowing down" of the bar is an artifact of 
unrealistically cold disk and grossly insufficient force resolution. 



the low and high resolution runs of the model A2 supports 
this picture. 

Numerical effects are clearly important for the structure 
of the system. Still, they can no t explain all the differences 
betwe en our results and those of lAthanassoula fc Misirioti^ 
l)2002h . The bar in the low resolution run of A2 model is 
longer than in the high resolution run. Still, it is shorte r 
than what was found bv lAthanassoula fc Misiriotij 1)200211 . 
The remaining differences are likely due to the differences 
in the overall structure of the system. Our halo was sig- 
ni ficantly more extended (250 kpc a s compared with 50 kpc 
in lAthanassoula fc MisiriotisI (|2002h 'l. Thus, the dark matter 
particles have larger velocities. Our disk is also slightly "hot- 
ter". Each effect points in the same direction: systems stud- 
ied by^;hanassoula & Misiriotis ( 2002) are likely more sus- 
ceptible to resonances. The effec ts accumulate over billions 
of yea r s producing lon ger bars in lAthanassoula fc Misiriotij 
l)2002h . lAthanassoulal jJoOS) presents a similar and more ex- 
tensive discussion which conclude that a hotter system is less 
resonant and has a slower decline of the bar pattern speed 
and a lower rate of bar growth. In particular lAthanassoulal 
l|2003tl shows that bars in disks with a Toomre parameter 
Q = .1 slow down more thant four times faster than in 
a disk with Q = 2.2, this result is in agreement with our 
discussion. We can conclude that the structure and kine- 
matics of halo and disk have a central role in bar evolution. 
Differences within our sim ulations and previous works like 
lAthanassoula fc Misirioti^ ( i2002) ) are influenced not only 



Figure 17. Effects of resolution on bar formation in model A2. 
Left panels show stellar particles in the simulation, which reached 
20 pc resolution. Right panels show the same model simulated 
with the resolution limited to 350 pc. The top panels are for 
2 Gyrs after initial moment; the bottom panels are for 4.5 Gyrs. 
The high resolution run has a very short (1 kpc) bar for 3 Gyrs, 
while the low resolution run has remarkably strong and long bar 
from the very beginning. Axes are labeled in units of kpc. This 
figure illustrates that numerical resolution can significantly affect 
structure of a bar. Lack of resolution results in longer and stronger 
bars. 



by resolution but also by the different structure of the mod- 
els. 

ft is still a matter of debate which initial setup 
is more realistic. Our halos are consistent with predic- 
tions of cosmological models. They extend to the virial 
radius; t hey have realistic profiles and r otation curves. 
Halos in lAthanassoula fc Misirioti^ i2002l) produce rea- 
sonable rotation curves although they are not consistent 
with cosmological models, fn an y case, we do not confirm 
lAthanassoula fc Misirioti^ (|22q3) conclusions that bars in 
models with centrally concentrated dark matter are neces- 
sarily too long (3-4 disk initial exponential lengths). Their 
results are valid only for the particular models they stud- 
ied. They are not valid as a general statement that re- 
alistic dark matter models should have excessively long 
bars. While we have men t ioned some disagreements with 
lAthanass oula fc Misirioti^ (|20o3), it should bo emphasized, 
that we find many similar r esults at least o n the q ualita - 
tive level. In agreement w ith Dcbattis ta fc SellwoodI (l200Cl) 
and lAthanassoula fc Misiriotisn (120021) we find that a mas- 
sive dark matter halo does not prevent formation of a strong 
bar. We also find that bars have a tendency to grow with 
time and do not dissolve spont aneously. We agree with 
lAthanassoula fc Misirioti^ i2002l) that formation of a bar 
leads to a baryon-dominated central region even in cases 
when the disk is initially submaximal. We also find that for- 
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mation and evolution of bars does not lead to a decline in 
the density of the dark matter in the central region. 

It is an interesting question how a bar can possibly 
not slow down when it rotates inside a dense dark mat- 
ter halo. Even in model Ai the disk lost 6 percent of its 
angular momentum after 5 Gyrs of evolution. One should 
expect that most of the loss comes from the bar. Because 
the bar has only 40% of the stellar mass, one expects that 
the bar should lose around 15% of its angular momentum. 
This is still marginally consistent with the behavior of the 
pattern speed in Figure ITI11 Still, what if Qp is really con- 
stant? Can a bar rotate with the constant speed and lose 
its angula r moment um? If we believe the arguments pre- 
sented bv IWeinberd (jlOSfil ). the answer is "No". Indeed, if 
a solid rotating bar, often used in these kind of estimates 
JWeinberdlioia iLittle fc CarlberdllQQlt IWeinberg fc Katj 
120011 ^. loses its angular momentum, it must slow down. Yet, 
real bars may behave differently because they are not solid 
objects. Bars are made of particles, which stream through 
them creating a wave pattern, which we call a "bar" . When 
the bar loses angular momentum, it is individual particles 
which lose it. The particles can move to orbits with smaller 
radii where orbital frequency is larger. It is reasonable to 
assume that the bar speed is proportional to the orbital fre- 
quency of particles. Because the particles rotate with larger 
frequency, the bar can speed up when it loses the angular 
momentum. Real bars have another feature, which compli- 
cates things even more. Typically, they have a tendency to 
grow. For example, the bar in model Ai has increased its 
length from ~ 4 — 5 kpc at 3 Gyrs to « 6 — 6.5 kpc at 5 Gyrs. 
This 30% increase in the bar length is due to particles, that 
came from large radii. Because the rotation curve is nearly 
flat, particles at larger radii have smaller orbital frequen- 
cies. When they "join" the bar, it will have a tendency to 
slow down. This happens if the particles do not significantly 
change trajectories. If the particles lose their angular mo- 
mentum and decrease average radius, they may cause the 
bar to speed up. 

The outcome of all these competing processes is difficult 
to predict. We find that slow growing bars rotate with almost 
constant speed (e.g., bars in models Ai, B, and the bar in 
model A2 at t < 3 Gyrs). If a bar grows fast, it significantly 
decreases its pattern speed. This happened in model A2 at 
t = 3 — 5 Gyrs. At 3 Gyr the bar was very short (~ 1.5 kpc) 
and weak. At earlier moments of time it was growing very 
slowly and was rotating with high constant speed for almost 
3 Gyrs. At t ~ 3 Gyrs the bar starts to grow very fast. By 
the end of evolution (6 Gyrs) it was 4.5 - 5 kpc long and it 
slowed down by a factor of two. 

Do bars in our models rotate as fast as the ob- 
served bars? In our models bars have TZ — 1.2 — 1.7. 
This should be compared with the observational results. 
Unfortunately the number of galaxies, for which the di- 
rect method of .Tremaine fc Weinberg 1^3) used, is 
small and the errors are s till la rge. For NGC 936 (Sa 
type) iMerrifield fc Kuiikeiil lll995l) find 7^ = 1.4 ± 0.3. 
iGerssen et alJ (|l99flD find ^ = l.lSl^-g^ fo r SO ga laxy NGC 
4596. For five SO galaxies lAguerri et aP 1)20031) find that 
"TZ is consistent with being in the range from 1.0 to 1.4". 
TaJien at the face value, these results are not consistent with 
most of our bars, which typically have TZ « 1.5. Yet, on 
pure statistical grounds our bars are still acceptable because 



the errors of the quoted observational results are only la's. 
For example, at the 2a level, the five galaxies studied by 
lAguerri et all l)2003L TableS) give the following upper lim- 
its: TZ < 1.4 (ESQ 139-G009), < 2.8 (IC 874), < 1.6 (NGC 
1308), < 1.9 (NGC 1440), and < 2.3 (NGC 3412). Our bars 
are well within the limits. 

There is another alarming issue: the overall structure 
of early type - SO and Sa - spirals studied in the observa- 
tions is very different from that of the numerical models. The 
later aim at Sb and Sc types, not SO. The surface density in 
simulations is well described by the double exponential law 
(Figure 15) ; the stellar velocity dispersion is relatively small 
(Figures 2- 3). The situation is v ery different for the galaxies 
studied bv lAguerri et al.l i2003l) . Surface brightness profiles 
are not even close to double exponentials. Every galaxy re- 
quires four comp onents: a central bul ge, a bar, a lens, and 
a disk (Figure 7. lAguerri et alJl2003h . For example, NGC 
3412, which has circular velocity 205kms~^ (close to our 
models) , is dominated in the central ~ 1 kpc region by a 
bulge, which does not have any counterpart in the simu- 
lations. A massive lens component is important at around 
3 kpc. Again it does not exist in our models. Line-of-sight 
velocities are also significantly larger than what we have in 
the models. While it is clear that we are dealing with very 
different objec ts, one may only speculate h ow these differ- 
ences affect TZ. ICombes fc ElmegreenI il993h find that early 
type spirals should have faster rotating bars. If that is con- 
firmed for models with live bulges and halos, disagreements 
with observations m ay become sma l ler. In the same direc- 
tion point results of lAthanassoulal (|20o3), who finds that 
bars in hotter disks rotate faster. 

In order to resolve the issues, numerical models should 
be improved to mimic SO galaxies. Models likely should in- 
clude a small initial bulge and should have disks with larger 
random velocities. Observations should also be done for Sb 
and Sc galaxies. Before this is done, it is difficult to assess 
the situation. As it stands now, differences with observations 
are not large and models cannot be rejected. 



8 CONCLUSIONS 

We study formation of bars in galactic models with equilib- 
rium exponential disks rotating inside extended dark mat- 
ter halos predicted by modern cosmological models. The 
amount of the dark matter in the central 10 kpc region is 
substantial - it is larger than the disk mass. Our models have 
parameters, which approximately match those of the Milky 
Way galaxy: the virial mass is (1 — 2) x 10^^ Mg , the total 
mass of the disk and the bulge is (4 — 6) x 10^^ Mq . Dark 
matter halos in our models extend to the virial radius of 
200-300 kpcand have cuspy inner profiles. We use millions 
of particles to simulate the models. 

Extended strong bars form in all our models in spite 
of the fact that the random velocities are much larger than 
the kinetic energy of disk rotation. Our results are at odds 
with the Ostriker-Peebles criterion of bar stability. We find 
that the dynamical friction between the dark matter and 
the bar results in some transfer of angular momentum to 
the dark matter halo, but the effect is much smaller then 
previously found in low resolution simulations and is incom- 
patible with early analytical estimates. The mass and the 
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force resolutions arc crucial for the dynamics of bars. In sim- 
ulations with the low resolution of 300-500 pc we find that 
the bar slows down and the angular momentum is lost rel- 
atively fast. In simulations with millions of particles, which 
reach the resolution of 20-40 pc, the pattern speed may not 
change over billions of years and the stellar component loses 
very little (5 - 10%) of its total angular momentum. 

Bars in our models are fast rotators with the ratio of 
the corotation radius to the bar major semi-axis being in 
the range TZ = 1.2 — 1.7, which is marginally compatible 
with observational results. The initial random velocities of 
the disk in the central fa 1 — 2 kpc region can substantially 
affect the structure of bars. In our simulations "hot" disks 
develop a longer bar, which rotates at an almost constant 
pattern speed. Colder models generate an initially shorter 
bar, which grows and slows down. 

In contrast to many previous simulations, bars in our 
models are relatively short. As in many observed cases, the 
major bar semi-axis in the models is about an exponential 
length of the disk. 

The transfer of the angular momentum between inner 
and outer parts of the disk plays a very important role in 
the secular evolution of disk. The main effect is the angular 
momentum transfer from the central part of the bar to the 
outer disk. This transfer leads to a substantial increase of 
the stellar mass and to a decrease of the dark matter-to- 
baryons ratio in the centre of the galaxy. For the Milky Way- 
scale models the central 2 kpc is strongly dominated by the 
baryonic component after 1-2 Gyrs since the onset of the 
instability. At intermediate 3 -10 kpc scales the disk is sub- 
dominant: spherically average density of the disk is 2-3 times 
smaller than the density of the dark matter. This makes 
an interesting twist for the never ending debate between 
supporters and opponents of the maximum disks. We predict 
that barred galaxies (or galaxies, which were barred) in the 
central region ( ~ 1/2 of the disk exponential length) should 
have the maximum disk (maximum bulge or bar, to bo more 
precise), but the outer part is the sub- maximum disk. 

Another effect of the bar formation is the increase by a 
factor of w 1.2 — 1.5 of the exponential length of the disk. 
This increase may change theoretical predictions for the disk 
lengths in hierarchical models. 

We find that the surface density profile of an evolved 
system is well approximated by a double exponential law. 
The 1/4 law for the bulge gives a worse fit, but it is not ex- 
cluded. Only in extreme cases (very late stages of evolution, 
very strong bars formed in low resolution simulations) we 
see a tendency for a bar to produce a flat part of the surface 
density in the interface of the bar and the disk. 

To summarize, we use more realistic models than in 
most of previous simulations and find that the models with 
substantial amount of the dark matter produce bar structure 
in striking agreement with observational results. 
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